
# Load packages
  library(tidyverse)

# Prepare data

  # Load Data-set
  dt <- readRDS("ObsStudy_00_data_Sweden_2018.RDS")
  
  # Vote 
  df <- dt %>% select(contains("ptv")) 
  dt$vote <- factor(names(df)[max.col(df)],
                    levels = c("ptv_M","ptv_SAP","ptv_C","ptv_V","ptv_SD","ptv_KD","ptv_MP","ptv_L"),
                    labels = c("M","SAP","C","V","SD","KD","MP","L"))
  vote_SAP <- dt$ptv_SAP
  vote_M <- dt$ptv_M
  
  # Coalition Evaluation
  Z <- dt %>% select(contains("rat_coal")) %>% mutate_all(funs(scales::rescale(.,to = c(0, 1), from = c(1,7)))) 
  rat_coal_SAP <- Z %>% select(rat_coal_SAPMP,rat_coal_SAPMPLC)
  rat_coal_M <- Z %>% select(rat_coal_ALLIANCE,rat_coal_MSD) 
  
  # Coalition Liklihood
  gamma <- dt %>% select(contains("coallik_")) 
  coal_lik_SAP <- gamma %>% select(coallik_SAP_MP,coallik_SAP_MPLC) %>%
    mutate(sum_exp = exp(coallik_SAP_MP) + exp(coallik_SAP_MPLC),
           coallik_SAP_MP = exp(coallik_SAP_MP)/sum_exp,
           coallik_SAP_MPLC = exp(coallik_SAP_MPLC)/sum_exp) %>%
    select(-sum_exp)
  
  coal_lik_M <- gamma %>% select(coallik_M_ALLIANCE,coallik_M_SD) %>%
    mutate(sum_exp = exp(coallik_M_ALLIANCE) + exp(coallik_M_SD),
           coallik_M_ALLIANCE = exp(coallik_M_ALLIANCE) / sum_exp,
           coallik_M_SD = exp(coallik_M_SD) / sum_exp) %>%
    select(-sum_exp)

# Mean-Variance Model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  coal_lik_M <- as.matrix(coal_lik_M)
  rat_coal_M <- as.matrix(rat_coal_M)
  coal_lik_SAP <- as.matrix(coal_lik_SAP)
  rat_coal_SAP <- as.matrix(rat_coal_SAP)
  
  M_M <- sapply(1:nrow(coal_lik_M), function(i) EV(coal_lik_M[i,], rat_coal_M[i,]))
  V_M <- sapply(1:nrow(coal_lik_M), function(i) Vari(coal_lik_M[i,], rat_coal_M[i,]))

  M_SAP <- sapply(1:nrow(coal_lik_SAP), function(i) EV(coal_lik_SAP[i,], rat_coal_SAP[i,]))
  V_SAP <- sapply(1:nrow(coal_lik_SAP), function(i) Vari(coal_lik_SAP[i,], rat_coal_SAP[i,]))

  # Data Frame
  d_s <- cbind(dt, 
             M_M,
             V_M,
             M_SAP,
             V_SAP
             )
  
  d_s$ptv_M <- scales::rescale(d_s$ptv_M, to = c(0, 1), from = c(1,7))
  d_s$ptv_SAP <- scales::rescale(d_s$ptv_SAP, to = c(0, 1), from = c(1,7))
  d_s$V_M <- scales::rescale(d_s$V_M, to = c(0, 1), from = c(0,0.25))
  d_s$V_SAP <- scales::rescale(d_s$V_SAP, to = c(0, 1), from = c(0,0.25))

# Saving model results

  # Moderates
  m1_s_tidy <- tidy(m1_s <- lm(ptv_M ~ M_M + rat_M + sex + as.factor(edu) + as.factor(age), d_s))
  m2_s_tidy <- tidy(m2_s <- lm(ptv_M ~ M_M + V_M + rat_M + sex + as.factor(edu) + as.factor(age), d_s))

  SE_2018_OS_M_Variance_Estimate <- m2_s_tidy %>% filter(term=="V_M") %>% select("estimate")
  SE_2018_OS_M_Variance_SE <- m2_s_tidy %>% filter(term=="V_M") %>% select("std.error")
  SE_2018_OS_M_Mean_Estimate <- m2_s_tidy %>% filter(term=="M_M") %>% select("estimate")
  SE_2018_OS_M_Mean_SE <- m2_s_tidy %>% filter(term=="M_M") %>% select("std.error")
  
  # Social Democrats
  m5_s_tidy <- tidy(m5_s <- lm(ptv_SAP ~ M_SAP + rat_SAP + sex + as.factor(edu) + as.factor(age), d_s))
  m6_s_tidy <- tidy(m6_s <- lm(ptv_SAP ~ M_SAP + V_SAP + rat_SAP + sex + as.factor(edu) + as.factor(age), d_s))

  SE_2018_OS_SAP_Variance_Estimate <- m6_s_tidy %>% filter(term=="V_SAP") %>% select("estimate")
  SE_2018_OS_SAP_Variance_SE <- m6_s_tidy %>% filter(term=="V_SAP") %>% select("std.error")
  SE_2018_OS_SAP_Mean_Estimate <- m6_s_tidy %>% filter(term=="M_SAP") %>% select("estimate")
  SE_2018_OS_SAP_Mean_SE <- m6_s_tidy %>% filter(term=="M_SAP") %>% select("std.error")

